
%% Preliminary
ResultName  =   'Result_Baseline';
Folder = ['TableGraphs\\IRF_Decomposition\\'];

load([ResultName,'.mat']);
mkdir(Folder);

%% IRF
% StdVec      =   [1 1 1 -1 -1 1]'*0.01/4;
% T           =   50;
% IRF         =   IRF_1order(MODEL,solution,StdVec,[],T);
T           =   50;
%% Decompose the IRF
VarList     =   {'Dist','Pir','T','Profit','TAU','ir','w'};
ShockList   =   {'Eps_Eta_M','Eps_Eta_Z'};
InputDev    =   struct();
for ii=1:length(ShockList)
    SH          =   ShockList{ii};
    for jj=1:length(VarList)
        VV          =   VarList{jj};
        if isfield(IRF.(SH),VV)
            InputDev.(SH).(VV)  =   IRF.(SH).(VV);
        elseif isfield(IRF.(SH),VV(1:end-1))
            Temp                =   zeros(size(IRF.(SH).(VV(1:end-1))));
            Temp(:,1:end-1)     =   IRF.(SH).(VV(1:end-1))(:,2:end);
            InputDev.(SH).(VV)  =   Temp;
        end
    end
end

DecDev      =   struct();
for ii=1:length(ShockList)
    SH          =   ShockList{ii};
    DecDev.(SH) =   SubFun_UnitIRFDec(PP,SS,MODEL,EquJac,VarList,T,InputDev.(SH));
end

%% Group-wise Decomposition
% Setup
VarGroup    =   struct('Interest',{{'ir'}},'Wage',{{'w'}},'Inflation',{{'Pir'}},...
                       'TAU',{{'TAU'}},'T',{{'T'}},'Profit',{{'Profit'}}, ...
                       'Distribution',{{'Dist'}});
VarGroupLabel=  struct('Interest','Interest Rate','Wage','Wage','Inflation','Inflation', ...
                       'TAU','Tax','T','Transfer','Profit','Dividend',...
                       'Distribution','Distribution'...
                       );
GroupList   =   fields(VarGroup);
GroupLabelList= cell(length(GroupList),1);
for ii=1:length(GroupList)
    GroupLabelList{ii}  =   VarGroupLabel.(GroupList{ii});
end
GroupDecDev =   struct();
for ii=1:length(ShockList)
    SH          =   ShockList{ii};
    TempList    =   fields(DecDev.(SH));
    for jj=1:length(TempList)
        VV          =   TempList{jj};
        GroupDecDev.(SH).(VV) ...
                    =   zeros(MODEL.VAR.(VV).Dim,T,length(GroupList));
        for gg=1:length(GroupList)
            GG          =   GroupList{gg};
            GGTempList  =   VarGroup.(GG);
            for kk=1:length(GGTempList)
                GroupDecDev.(SH).(VV)(:,:,gg)    =   GroupDecDev.(SH).(VV)(:,:,gg)+DecDev.(SH).(VV).(GGTempList{kk});
            end
        end
    end
end

%% Collect the Results
TempVarList =   {'C','AggStat'};
IrfAgg      =   struct();
IrfDec      =   struct();
IrfDecTable =   struct();
Table_Rows  =   GroupList;
Table_Cols  =   {'C';...
                 'poor';'rich'};
for ii=1:length(ShockList)
    SH          =   ShockList{ii};
    % Peak IRF of Aggregate C
    [~,AggInd]  =   max(abs(IRF.(SH).C));
    % First 8 Quarters
%     AggInd      =   (1:8);
    % Normalizer
    Normalizer  =   abs(sum(IRF.(SH).C(:,AggInd)));
    for jj=1:length(TempVarList)
        vv          =   TempVarList{jj};
        TempSize    =   size(GroupDecDev.(SH).(vv));
        IrfAgg.(SH).(vv) ...
                    =   reshape(sum(GroupDecDev.(SH).(vv)(:,AggInd,:),2),TempSize([1,3])) ...
                        /Normalizer;
        IrfDec.(SH).(vv) ...
                    =   IrfAgg.(SH).(vv)./sum(IrfAgg.(SH).(vv),2);
    end
    %----------------------------------------------------------------------
    % Tables
    %----------------------------------------------------------------------
    % Aggregated
    Table_Data  =   [IrfAgg.(SH).C;IrfAgg.(SH).AggStat]';
    IrfAggTable.(SH) ...
                =   array2table(Table_Data,'VariableNames',Table_Cols,'RowNames',Table_Rows);
    IrfAggTable.(SH).Diff_poor_rich ...
                =   IrfAggTable.(SH).poor-IrfAggTable.(SH).rich;
    % Pct
    IrfDecTable.(SH) ...
                =   IrfAggTable.(SH);
    for jj=1:length(IrfAggTable.(SH).Properties.VariableNames)
        vv          =   IrfAggTable.(SH).Properties.VariableNames{jj};
        IrfDecTable.(SH).(vv) ...
                    =   IrfDecTable.(SH).(vv)/sum(IrfDecTable.(SH).(vv));
    end
    % Print Tables
     for jj=1:length(IrfAggTable.(SH).Properties.VariableNames)
        vv          =   IrfAggTable.(SH).Properties.VariableNames{jj};
        IrfAggTable.(SH).(vv)   =   round(IrfAggTable.(SH).(vv),3);
        IrfDecTable.(SH).(vv)   =   round(IrfDecTable.(SH).(vv),3);
    end
    writetable(IrfAggTable.(SH),[Folder,'IrfAgg.xlsx'],'WriteRowNames',true,'Sheet',SH);
    writetable(IrfDecTable.(SH),[Folder,'IrfDec.xlsx'],'WriteRowNames',true,'Sheet',SH);
end


%% Visual Decomposition by Groups
N_col       =   1;
N_row       =   1;
fig_size    =   [N_col,N_row].*[0.75,1]*1/2; 
fig_gap     =   [0.05,0.08];
fig_VMargin =   [0.05,0.03];
fig_HMargin =   [0.10,0.02]; 

fig_SubPlot =   @(ii_vv)subtightplot(N_row,N_col,ii_vv,...
                                     fig_gap,fig_VMargin,fig_HMargin);
AxFontSize  =   8;
TickFmt     =   '%.2f';
XLabel      =   @(Label)xlabel(Label,'fontsize',8,'interpreter','latex');
YLabel      =   @(Label)ylabel(Label,'fontsize',8,'interpreter','latex');
LEG         =   @(PH,LabelList)legend(PH,LabelList,'fontsize',8,'interpreter','latex','location','northeast');


for ii=1:length(ShockList)
    SH = ShockList{ii};
    
    YNames = cell(length(IrfAggTable.(SH).Row),1);
    for jj=1:length(YNames); YNames{jj} = VarGroupLabel.(IrfAggTable.(SH).Row{jj}); end
    % Financial Integration
    % by Net Wealth
    TempData = [IrfAggTable.(SH).poor,IrfAggTable.(SH).rich,IrfAggTable.(SH).Diff_poor_rich];
    XNames = {'Low Net Asset','High Net Asset','Difference'};

    Fig = FigureSetup([SH,', Net Asset'],fig_size);
    ax = fig_SubPlot(1);
    [~,ax] = Plot_StackedBars(TempData,YNames,XNames,[],ax);
    
    print('-depsc','-r1000',Fig,[Folder,'byNA_',SH]);
    close(Fig);
end

